clear all
set mem 1000
cd "D:\data\origin\COVID19\"
use germanydata


//Fig. 3

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick)) 
	graph save "germanyresult/FIGURE3_1.gph", replace
	restore
	
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths  east , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE3_2.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio  east , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium)) ylab(,labsize(vsmall)) ylabel(-2(1)2) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE3_3.gph", replace
	restore	
	
//Fig. 5
	local typelist "border"	
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(small)) xtitle("") title("Districts on the Border",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/cases`type'.gph", replace
	restore
		
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(small))   xtitle("") title("Districts on the Border",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deaths`type'.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(small))   xtitle("") title("Districts on the Border",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deathr`type'.gph", replace
	restore			
	}
	
	local typelist "distance50"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(small)) xtitle("") title("Districts within 50 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/cases`type'.gph", replace
	restore
	
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(small))   xtitle("") title("Districts within 50 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deaths`type'.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(small))   xtitle("") title("Districts within 50 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deathr`type'.gph", replace
	restore
	}

	local typelist "distance100"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(small)) xtitle("") title("Districts within 100 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/cases`type'.gph", replace
	restore
	
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(small))   xtitle("") title("Districts within 100 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deaths`type'.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(small))   xtitle("") title("Districts within 100 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deathr`type'.gph", replace
	restore
	}

	local typelist "distance150"
		
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(small)) xtitle("") title("Districts within 150 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/cases`type'.gph", replace
	restore
	
	
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(small))   xtitle("") title("Districts within 150 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall)) ylabel(-1(0.5)1) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deaths`type'.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio east if `type'==1, robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen id2=_n-275
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(small))   xtitle("") title("Districts within 150 km",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medsmall)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/deathr`type'.gph", replace
	restore
	}
	
	gr combine "germanyresult/casesborder" "germanyresult/casesdistance50" "germanyresult/casesdistance100" "germanyresult/casesdistance150" , saving("germanyresult/FIGURE5_1",replace) 
	gr combine "germanyresult/deathsborder" "germanyresult/deathsdistance50" "germanyresult/deathsdistance100" "germanyresult/deathsdistance150" , ycommon saving("germanyresult/FIGURE5_2",replace) 
	gr combine "germanyresult/deathrborder" "germanyresult/deathrdistance50" "germanyresult/deathrdistance100" "germanyresult/deathrdistance150" , ycommon saving("germanyresult/FIGURE5_3",replace) 


	
	
//Fig. 6

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  vaccination  east , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20210101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Vaccination)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2021"  105"Apr 2021" 196"Jul 2021" 288"Oct 2021" 349 "Dec 2021"  , labsize(medium)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE6_1.gph", replace
	restore

	local typelist "border"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  vaccination  east if `type'==1 , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20210101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Vaccination)",size(small)) xtitle("") title("Districts on the Border",size(medlarge)) xlabel(15 "Jan 2021"  105"Apr 2021" 196"Jul 2021" 288"Oct 2021" 349 "Dec 2021" , labsize(small)) ylab(,labsize(vsmall)) ylabel(-30(10)10)  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/vacc`type'.gph", replace
	restore
	}
	
	local typelist "distance50"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  vaccination  east if `type'==1 , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20210101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Vaccination)",size(small)) xtitle("") title("Districts within 50km",size(medlarge)) xlabel(15 "Jan 2021"  105"Apr 2021" 196"Jul 2021" 288"Oct 2021" 349 "Dec 2021"  , labsize(small)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/vacc`type'.gph", replace
	restore
	}
	
	local typelist "distance100"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  vaccination  east if `type'==1 , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20210101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Vaccination)",size(small)) xtitle("") title("Districts within 100km",size(medlarge)) xlabel(15 "Jan 2021"  105"Apr 2021" 196"Jul 2021" 288"Oct 2021" 349 "Dec 2021"  , labsize(small)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/vacc`type'.gph", replace
	restore
	}

	local typelist "distance150"
	
	foreach type in `typelist' {
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  vaccination  east if `type'==1 , robust
	rename _eq2_stat_1 number
	replace number=number-1
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20210101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Vaccination)",size(small)) xtitle("") title("Districts within 150km",size(medlarge)) xlabel(15 "Jan 2021"  105"Apr 2021" 196"Jul 2021" 288"Oct 2021" 349 "Dec 2021"  , labsize(small)) ylab(,labsize(vsmall))  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/vacc`type'.gph", replace
	restore	
	}
	
	gr combine "germanyresult/vaccborder" "germanyresult/vaccdistance50" "germanyresult/vaccdistance100" "germanyresult/vaccdistance150" , ycommon saving("germanyresult/FIGURE6_2",replace) 

//Fig. 7

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  confirmedcases east vaccination land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police relativeday relativeday_2 , robust
	rename _eq2_stat_1 number
	replace number=number-20
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Confirmed Cases)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium)) ylab(,labsize(vsmall)) ylabel(-12(3)0)  graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE7_1.gph", replace
	restore		
			
	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deaths east vaccination land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police relativeday relativeday_2, robust
	rename _eq2_stat_1 number
	replace number=number-20
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Cases)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium))  ylab(,labsize(vsmall)) ylabel(-10(2)0) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE7_2.gph", replace
	restore

	preserve
	sort newdate
	statsby _b _se e(N), by(newdate) clear:reg  deathratio east vaccination land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police relativeday relativeday_2, robust
	rename _eq2_stat_1 number
	replace number=number-20
	gen lower=_b_east-_se_east*invttail(number,0.05)
	gen upper=_b_east+_se_east*invttail(number,0.05)
	drop if newdate<20200101 | newdate>20211231
	sort newdate
	gen id=_n
	gen plotgroup= "2"
	save temp1,replace
	replace lower=_b_east-_se_east*invttail(number,0.025)
	replace upper=_b_east+_se_east*invttail(number,0.025)
	replace plotgroup= "1"
	save temp2,replace
	append using temp1
	eclplot _b_east lower upper id, estopts( msize(vsmall)) ytitle("Coefficient (Regression of Death Ratio)",size(medlarge)) xtitle("") title("",size(medlarge)) xlabel(15 "Jan 2020"  197"Jul 2020" 381"Jan 2021" 562 "Jul 2021" 715 "Dec 2021" , labsize(medium)) ylab(,labsize(vsmall))ylabel(-60(30)60) graphregion(color(white)) supby(plotgroup) ciopts1(blcolor(eltblue)) ciopts2(blcolor("78 155 235")) legend(off) yline(0, lcolor(red) lwidth(medthick))
	graph save "germanyresult/FIGURE7_3.gph", replace
	restore


//Table 3
	//to calculate mean and median
	local datelist "20200415 20201015 20210415 20211015"
	foreach dat in `datelist' {
	eststo clear
	eststo: qui estpost tabstat confirmedcases  confirmedcases_unlogged deaths deaths_unlogged deathratio vaccination if newdate==`dat'  , by (east) statistics(n mean p50) columns(statistics)
	esttab est* using germanyresult/TABLE3_1.csv, append  cells("count(lab(N)) mean(fmt(2)) p50(fmt(2))") 
	}
	eststo clear
	eststo: qui estpost tabstat land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police if newdate==20200415  , by (east) statistics(n mean p50) columns(statistics)
	esttab est* using germanyresult/TABLE3_2.csv, append  cells("count(lab(N)) mean(fmt(2)) p50(fmt(2))") 

	//to calculate t-value and Wilcoxon rank-sum z-value
	local datelist "20200415 20201015 20210415 20211015"
	local tlist "confirmedcases  confirmedcases_unlogged deaths deaths_unlogged deathratio vaccination"
	foreach dat in `datelist' {
	foreach var in `tlist' {
	ttest `var' if newdate==`dat'  , by (east) 
	ranksum `var' if newdate==`dat', by (east) 
	}
	}

	local tlist "land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police"
	foreach var in `tlist' {
	ttest `var' if newdate==20200415  , by (east) 
	ranksum `var' if newdate==20200415, by (east) 
	}	
	
	
//Table 4
	local ylist "confirmedcases deaths deathratio"
	local datelist "20200415 20201015 20210415 20211015"
	local cont "vaccination land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police relativeday relativeday_2"
	local numb "daily"
	
	foreach y in `ylist' {
	foreach dat in `datelist' {
	qui:reg `y' east if newdate==`dat', robust
	outreg2 using germanyresult/TABLE4_`y'.xls,append ctitle(`dat') sortvar(east `cont') bd(3) td(2) rd(3) tstat adjr2 excel
	qui:reg `y' east `cont'  if newdate==`dat', robust
	outreg2 using germanyresult/TABLE4_`y'.xls,append ctitle(`dat')  sortvar(east `cont') bd(3) td(2) rd(3) tstat adjr2 excel
	}
	}
	
	
//SUPPLEMENTAL TABLE 4

	local datelist "20200415 20201015 20210415 20211015"
	foreach dat in `datelist' {
	eststo clear
	eststo: qui estpost tabstat confirmedcases  confirmedcases_unlogged deaths deaths_unlogged deathratio vaccination relativeday relativeday_2 if newdate==`dat'  , statistics(n mean sd p10 p25 p50 p75 p90) columns(statistics)
	esttab est* using germanyresult/TABLES4_1.csv, append  cells("count(lab(N)) mean(fmt(2)) sd(fmt(2)) p10(fmt(2)) p25(fmt(2)) p50(fmt(2)) p75(fmt(2)) p90(fmt(2))") title(`dat')
	}
	eststo clear
	eststo: qui estpost tabstat  east land_area population_age65 grp_growth export healthcare migration  serviceperson student educationexp cinema residentialbuilding passenger accommodation  rail police if newdate==20200415  , statistics(n mean sd p10 p25 p50 p75 p90) columns(statistics)
	esttab est* using germanyresult/TABLES4_2.csv, append  cells("count(lab(N)) mean(fmt(2)) sd(fmt(2)) p10(fmt(2)) p25(fmt(2)) p50(fmt(2)) p75(fmt(2)) p90(fmt(2))")  title(District Characteristics)
